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is  introduced  into  the  analysis  in  order  to  facilitate  the  computational  task.  Although 
the  turbulent  flow  problem  is  inherently  nonlinear,  the  method  of  analysis  described 
herein,  typically  employed  with  linear  electrical  networks,  is  found  to  converge  rapidly 
to  an  approximate  solution  without  user-specified  a  priori  guesses. 
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1.  Introduction 

An  interest  in  the  properties  of  incompressible  fluid  flow  in  pipes  cam  undoubtedly  be 
traced  as  far  bamk  as  the  era  of  the  Roman  empire.  A  brief  synopsis  of  the  historical 
development  of  fluid  mechanics  can  be  foimd  in  many  of  the  introductory  level  texts 
[Vennard  (1961)  for  example]. 

The  problem  of  accurately  solving  for  the  parameteiTS  in  incompressible  fluid  flow 
in  piping  networks  is  dependent  upon  the  single- pipe  model  for  fluid  flow.  Many 
semi-empirical  methods  for  characterizing  the  flow  properties  have  appeared  in  the 
literature  [Matthew  (1981),  Wuori  (1993),  and  Churchill  (1977)].  Since  the  early  work 
by  Feigenbaum  (1978),  see  for  example  [Hofstadter  (1981)],  on  chaos  in  simple  sys¬ 
tems,  considerable  progress  has  been  made  on  analytically  characterizing  turbulence 
in  fluids  [Landau  and  Lifshitz  (1987)].  Nevertheless,  many  engineering  problems  are 
solved  using  the  Stanton  diagram,  which  shows  Nikuradse’s  measured  data. 

For  engineering  problems  which  require  repetitive  numerical  evaluation,  an  empir¬ 
ical  relationship  facilitates  the  task  and  various  investigators  [Tsang  and  Kee  (1987), 
Round  (1980,  1985),  and  Colebrook  (1939)]  have  proposed  such  formulae.  These  for¬ 
mulae  have  the  drawback  that  the  transition  region  between  laminar  and  turbulent 
flow  is  inadequately  represented.  Only  the  fluid  flow  properties  in  the  laminar  region 
are  fairly  well  understood  [Vennard  (1961)]  using  methods  of  fluid  analysis. 

One  of  the  earliest  and  most  well  known  methods  for  solving  pipeline  networks 
is  attributed  to  Hardy  Cross  [Giles  (1962)  and  Lindeburg  (1982)].  The  Hardy  Cross 
scheme,  not  only  requires  an  a  priori  guess  on  the  initial  flows,  but  upon  comparison 
with  more  recently  reported  methods  [Carnahan  and  Christensen  (1972),  Bending 
and  Hutchinson  (1973),  Gay  and  Preece  (1975,  1977)],  it  is  found  to  be  less  efficient. 
Finally,  in  a  recent  Naval  Postgraduate  School  (NPS)  thesis  [Ellis  (1993)],  a  pipeline 
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method  formalism  based  on  the  Cholesky  method  lor  matrix  reduction  has  been 
investigated.  Although  the  mathematical  analysis  of  the  scheme  is  fairly  efficient,  the 
physical  modeling  did  not  encompass  conditions  for  other  than  laminar  flows. 

The  main  purpose  of  this  report  is  to  extend  the  range  of  flow  to  include  transition 
and  turbulent  flow  regimes.  As  this  is  primarily  a  feasibility  study,  various  facilitating 
assumptions  were  made.  First,  it  is  assumed  that  English  units  are  preferred  for 
input/output.  Second,  all  pipes  are  assumed  to  have  uniform  circular  cross-sections. 

2.  Physical  Modeling 

For  incompressible  fluids,  two  of  Bernoulli’s  rules  [Lindeburg  (1982)]  for  conservation 
of  mass  and  energy  are,  the  continuity  equation 


Q  =  A,V,  =  AaV, 


(1) 


and  Bernoulli’s  law 


^  +  (P2  -  Pi)  +  pp(^2  -  ^i)  +  pgh' 


(2) 


where  the  subscript  2/1  refers  to  positions  2/1  and  the  symbol  lists  given  in  Section 
10  defines  the  other  relevant  parameters.  The  head  loss,  as  defined  by  the  Darcy- 
Weisbach  equation,  is  given  by 


2gd 


(3) 


where  L  is  the  pipe  length,  d  is  the  hydraulic  diameter,  and  /  is  the  friction  factor. 
For  circular  cross-section  pipes  the  hydraulic  diantMter  and  the  physical  diameter  are 
identical.  Using  similitude  analysis  [Vennard  (1961)],  it  is  possible  to  show  that  the 
friction  factor  depends  only  on  two  dimensionless  parameters,  the  Remolds  number 

pVd 


Re  = 


(4) 


2 


where  /x  is  the  dynamic  viscosity,  and  the  relative  roughness 


c  = 


c 

d 


(5) 


where  e  is  the  roughness  of  the  pipe. 

A  sourceless  section  of  pipe  which  is  horizontal  (z2  =  zi),  and  with  constant 
cross-section  so  that  V3  =  Vi,  will,  from  eq  (2)  and  eq  (3),  then  have  a  pressure  drop 


Pi  -  P2  =  h  =  pgh*  = 


pfLV^ 

2d 


(6) 


or  after  use  of  eq  (1) 


(7) 


This  suggests  the  useful  electrical  analogy  with  Ohm’s  law  where  pressure  drop  cor¬ 
responds  to  voltage  and  fluid  flow  corresponds  to  current.  It  follows  that  consistent 
with  the  analogy,  the  effective  resistance  of  the  pipe  is  then  given  by 

rLpV 


R^f 


2dA 


(8) 


where  the  units  of  this  resistance  are  defined  by  the  ratio  of  pressure  to  flow. 

It  is  observed  that  the  problem  of  predicting  fluid  flow  in  pipes  now  reduces  to 
accurate  characterization  of  the  friction  factor.  For  the  laminar  flow  domain  the 
analytically  derivable  [Vennard  (1961)]  result 


/  = 


W 

Re 


(9) 


agrees  very  well  with  measurement  out  to  Re  <  2100  which  is  the  well  known  limit  in 
pipes,  ducts,  tubes,  and  channels  for  laminar  flow.  After  substitution  of  eq  (9)  into 
eq  (8)  it  is  found  that: 


32Lft 

<PA 


(10) 
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and  the  corresponding  relationship  for  head  loss  is 


hi  = 


(11) 


This  is  in  agreement  with  an  expression  cited  by  Bending  and  Hutchison  (1993).  The 
subscript  t  has  been  used  to  indicate  the  laminar  flow  domain.  Note  from  eq  (10) 
that,  as  expected,  the  pipe  resistance  is  independent  of  the  flow  rate  in  the  laminar 
regime. 

In  order  to  easily  characterize  the  degree  to  which  the  fluid  flow  is  turbulent,  the 
head  loss  can  be  expressed  as 

h  =  A/  X  0  (12) 


where  the  turbulence  factor,  denoted  0,  is  a  linear  gauge  of  the  deviation  of  the 
system  from  laminar  flow.  After  noting  that  eq  (7)  can  be  reexpressed  as 


(13) 


it  can  be  concluded  from  eq  (4)  that,  in  agreement  with  the  definition  of  the  Reynolds 
number,  the  turbulence  factors  satisfy 


^  ■'  64 


(14) 


As  seen  from  eq  (9)  V’  =  1  for  laminar  flow.  In  the  next  section  it  will  be  shown  that 
for  transition  and  turbulent  flow  V*  >  1- 

The  pipe  resistance  can  now  be  expressed  as 


R  =  Rtxil} 


(15) 


where  as  previously  noted,  Rt  does  not  depend  on  the  flow  rate. 


4 


3.  Friction  Factor  Estimation  for  Nonlaminar 
Fluid  Flow 

For  purposes  of  the  developmeut  here  it  is  necessary  to  introduce  two  distinct  non¬ 
laminar  regions.  First,  the  transition  region  satisfying 


Re„  >  Re  >  2100 


(16) 


where  Re^  is  usually  cited  typically  between  3500  and  4000  [Vennard  (1961)]  and  for 
the  turbulent  region 

Re  >  Re„  (17) 


One  of  the  most  commonly  employed  empirical  formulas  for  the  turbulent  flow  domain 
defined  by  eq  (17)  is  attributed  to  Colebrook  (1939) 


c  2.51 

3.7d  BevT 


(18) 


Although,  as  seen  in  eq  (18),  the  dependence  of  the  friction  factor  on  Reynolds  number 
is  implicit,  the  expression  will  converge  rapidly  upon  iterative  substitution. 

To  date  the  most  viable  predictors  for  the  transition  region,  see  eq  (16),  are 
empirical.  The  Stanton  diagram  [Vennard  (1961)]  shows  smoothly  varying  curves 
which  are  based  on  measurements  in  the  transition  region.  The  Moody  diagram, 
which  only  plots  expressions  for  eq  (9),  for  Re  <  2100,  and  eq  (18),  for  Re  >  Rco,  is 
not  useful  for  the  transition  region.  Bending  and  Hutchinson  (1973)  suggested,  in  a 
computer  analysis  of  piping  networks,  that  a  linear  interpolation  between  the  laminar 
and  the  turbulent  domain  be  made  to  provide  friction  data  in  the  transition  region. 
For  the  network  analysis  scheme  to  be  described  herein,  a  linear  approximation  for 
the  transition  region  was  foimd  to  be  insufficiently  accurate  upon  numerical  testing. 
In  particular,  the  iterative  program  with  the  linear  approximation  would  either  not 
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converge  or  converge  very  slowly.  The  problem  in  this  scheme  was  the  dramatic 
discontinuity  in  the  first  derivative  at  the  edges  of  the  transition  region.  A  greatly 
improved  modeling  scheme,  in  terms  of  convergence,  is  given  by 


/  = 


Re 

2100 


+ 


TT  (Re  -  2100) 


2  (Re^  -  2100)  J 


(19J 


where  /(Reo)  was  calculated  via  the  Colebrook  model  eq  (18).  Note  that,  like  the 
linear  approximation  model,  continuity  is  preserved  at  the  edges  of  the  transition 
region.  Because,  as  seen  from  eq  (19),  the  slope  at  the  transition  edges  is  zero,  the 
discontinuity  in  the  slope  is  much  less  dramatic  and  the  convergence  was  found  to  be 
much  more  rapid.  Figure  1  is  a  graphical  representation  for  these  curves  in  standard 
log- log  format  for  various  relative  roughness  factors. 

It  is  worth  noting  that  the  work  by  Churchill  (1977)  agreed  very  well  with  the 
scheme  just  described  for  Re  >  2100  and  provided  an  explicit  formula  for  the  friction 
factor,  /.  However,  for  Re  <  2100  (laminar  flow)  the  expression  deviated  significantly 
from  the  classical  expression  of  eq  (9).  Hence  it  was  not  used. 

The  argument  supporting  the  cldm  that  the  turbulence  factor  eq  (14)  is  generally 
greater  than  unity  is  loosely  based  on  the  following  observations.  First,  in  1913  Blasius 
proposed  an  empirical  relationship  for  smooth  pipe,  i.e.,  c  =  0  in  the  turbulent  flow 
regime  (Vennard  (1961)] 


/■BMoUl 


I»P« 


0.316 

Re»'^ 


3000  <  Re  <  100,000 


(20) 


Second,  the  construction  algorithm  for  predicting  /  in  the  transition  region,  eq  (19) 
requires  that 


(21) 


This  leads  to  the  inequality  that  is  established  from  a  comparison  of  eq  (20),  eq  (21), 
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and  Fig  1 


/rou^  pipe  ^  ftmooih  pipe  ^  fuauatu  —  (22) 

and  is  supported  by  analysis  and  measurement  [Vennard  (1961)].  The  inequality, 
>  1,  then  follows  directly  from  the  definition  of  the  turbulence  factor,  0,  gi^^en  by 
eq  (14). 

The  turbulence  factor  has  been  introduced  into  this  report  for  two  reasons.  Al¬ 
though  the  friction  factor  is  a  linear  indicator  for  the  system  response  of  the  fluid 
dynamics,  it  is  not  a  anear  indicator  for  the  deviation  of  the  system  from  laminar 
behavior.  Note  that  the  quantity  (Re  —  2100)  is  a  nonlintar  indicator  for  the  devia¬ 
tion  from  laminar  behavior.  The  purpose  of  this  work  has  been  to  extend  the  Ellis 
(1993)  effort  to  encompass  nonlaminar  regimes.  The  turbulence  factor  will  readily 
demonstrate  the  need  for  this  extension.  Second,  it  turns  out  that  the  laminar  pipe 
resistance,  iZ/,  can  be  computed  at  the  onset  of  the  problem  and  only  the  turbulence 
factors  need  to  be  updated  in  the  iterative  calculations  of  eq  (15)  rather  than  to  com¬ 
pute  eq  (8).  In  brief,  the  turbulence  factors  provide  some  computational  advantages 
which,  for  large  piping  networks,  would  be  significant. 

4.  Basic  Circuit  Modeling 

Following  the  electrical  analogy  for  an  arbitrary  pipe,  a  pressure  difference,  Apk, 
across  a  pressure  source,  Ap,k,  in  series  with  a  pipe  resistance,  rj^,  with  a  flow  qk, 


satisfies 


^Pk  =  qkTk  - 


consistent  with  the  convention  established  by  Fig  2.  The  corresponding  relation  in 
terms  of  a  current  source  is  given  by 

9k  =  9.k  +  VfcApfc  (24) 
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f-F ACTOR  MODEL 


^  o  ^  ra 

I  I 


(^OlDVi  NOIlDra£)0  TOO! 

Figure  1:  Friction  factor  model. 
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LOGIO(REYNOLD’S  NUMBER) 


+ 

:QAf.. 


Figure  2:  Pressure  source  circuit  model. 

where  yk  —  l/fk-  Figures  2  and  3  define  the  nomenclature  and  the  topology.  Equation 
(24)  can  be  compactly  represented  for  all  branches  or  pipes  in  matrix  format 

Q=sQ.+YAP  (25) 

where  it  is  seen  that  eq  (24)  and  eq  (25)  have  unresolved  unknowns  associated  with 
flow  rate  and  pressure.  The  condition  of  continuity  at  the  nodes 

J^qkj-0  ;  =  1,2, ...nr  (26) 

* 

where  nr  is  the  total  number  of  nodes  in  the  network,  provides  enough  infoimation  to 
resolve  the  flow  rates  and  pressure  drops.  A  formulation  for  this  process  is  presented 
in  the  next  section. 
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Figure  3:  Flow  source  circuit  model. 

5.  Mathematical  Formulation 

The  mathematical  formalism  is  presented  via  am  illustration  taken  from  the  Elhs 
(1993)  thesis.  The  piping  network  shown  in  Fig  4a  has  five  nodes  and  six  connecting 
branches.  It  is  also  seen  that  the  branch  between  nodes  1  and  5  contains  a  pump. 
This  pump  provides  a  source  for  fluid  flow  at  its  high  pressure  side  shown  at  node-1 
with  an  arrow.  The  low  pressme  side  is  the  corresponding  sink.  The  corresponding 
flow  graph  is  shown  in  Fig  4b.  Each  branch  now  has  an  arbitrarily  chosen  orientation 
arrow  which  defines  the  convention  for  positive  fluid  flow.  Note  that  consistent  with 
the  flow  direction  shown  for  branch  1,  the  source  shown  on  Fig  4a  would  be  assigned 
a  negative  numerical  vale.  The  first  step  in  the  mathematical  formalism  is  to  cre¬ 
ate  both  a  network  representation  for  the  fluid  pipeline  network  and  a  corresponding 
flow  graph.  The  flow  rates  in  the  b  branches,  in  this  case  6,  can  be  represented  as  a  bx  1 
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(a) 


Figure  4:  a)  Six  branch  pipeline  network;  b)  Associated  flow  graph  for  a)  [15]. 
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<  = 

1 

Branch  j  leaves  node  i 

(28a) 

4  = 

-1 

Branch  j  enters  node  i 

(28b) 

= 

0 

otherwise 

(28c) 

where  t  =  1, 2, . . . ,  nj>  nodes  and  j 

=  1,2, . . . ,  6  branches.  This  leads  to: 

1 

1  0  0  0  O' 

• 

0 

-11100 

C“  = 

0 

0-1010 

(29) 

0 

0  0-1-1  1 

-1 

0  0  0  0  -1 

where,  in  general,  C"  has  6  columns  and  nj  rows.  It  is  apparent  that  for  each  column, 


which  is  associated  with  a  node, 

k 

or  more  formally 


C“Q  =  0  (31) 

consistent  with  flow  rate  conservation  at  all  nodes  (see  eq  (26)).  In  order  to  reduce 
the  order  of  the  problem  by  1  it  is  useful  to  define  one  of  the  nodes  as  the  ground  or 
datum  reference.  In  principle  it  should  not  matter  which  node  is  taken  as  this  datum 
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node,  however,  in  the  coded  implementation  it  is  assumed  that  the  circuit  nodes 
are  numbered  such  that  the  datum  node  is  the  node  bearing  the  highest  numerical 
designator.  For  example,  as  seen  on  Figs  4a,  b,  the  node-5  is  the  datum  node.  Because 
the  assignment  of  nodes  is  independent  of  the  assignment  of  branches  and  is  arbitrary, 
there  is  no  sacrifice  in  the  generality  of  the  method.  Following  this  point,  the  node¬ 
branch  incidence  matrix,  C,  can  be  defined  according  to  the  rules  dictated  by  eq  (27) 
with  the  index  t  satisfies  the  condition  t  =  1, 2, . . .  ,n  where 

n  =  nr  —  1  (32) 

Thus,  for  this  example,  the  node-branch  incidence  matrix  is 

1  1  0  0  0  0 

0-1  1  1  00 

'"“0  0-1  0  10 

.0  0  0-1-11 

In  general  C  has  b  columns  and  n  rows.  It  follows  from  eq  (30)  that 

CQ  =  0  (34) 


(33) 


which  forces  a  condition  of  flow  rate  conservation  or  continuity  at  all  nodes. 

The  node  pressures,  defined  with  respect  to  the  reference  datum  node,  can  be 
represented  by  the  column  vector 


Pi 


Pn  J 


(35) 


which,  for  the  example  under  consideration,  the  node  pressure  vector  P  is  a  4  x  1 
column  vector  since  n  =  4.  The  node  branch  incidence  matrix  can  be  used  to  predict 
the  head  losses  from  the  rule 

H  =  C^’P  (36) 
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Given  the  matrix  Y„  and  the  source  vector  Q,,  any  number  of  matrix  inverting 
schemes  [Hamming  (1962)]  could  be  employed  to  evaluate  eq  (40).  As  suggested  in  the 
Ellis  (1993)  thesis  the  very  efficient  Cholesky  reduction  algorithm  is  applicable  here 
because  the  matrix  Yn  is,  not  only  symmetric,  but  positive  definite.  The  Cholesky 
reduction  into  upper  and  lower  triangular  matrices  permits  a  rapid  matrix  inversion 
using  Gaussian  elimination  as  described  in  Appendix  A. 

Once  the  node  pressures  given  by  eq  (40)  are  obtained  the  branch  head  losses  and 
flow  rates  are  easily  calculated  from  eq  (36)  and  eq  (39),  respectively.  Note  that,  once 
the  flow  rate  vector  is  obtained,  the  flow  rate  velocities  can  be  computed  according 
to  the  rule 


J 


(43) 


in  agreement  with  eq  (1).  The  corresponding  velocity  magnitudes,  (|ui |,  [ujI,  . . . ,  [ubI) 
can  then  be  applied  with  eq  (4)  to  calculate  the  Reynolds  numbers,  the  friction  factors, 
and  the  turbulence  factors  using  eq  (14).  The  pipe  resistances  are  calculated  from  eq 

(15)  and  the  corresponding  admittances,  yi,y2,. .  •  ,y6,  follow  from 

l/ri  0  0  •••  0 


y  = 


0 

0 


l/r. 


0 

0 

0 


(44) 


L  0  0  •  •  0  l/rt  J 

If  the  current  turbulence  factors  differ  significantly  from  the  previous  turbulence  fac¬ 
tors  according  to  the  rule 


MAXERR  >  max  |^*(previous)  -  V'*(current)|  :  {branches  ib  =  1, 6}  (45) 


the  process  is  repeated.  The  Fortran  input  parameter  MAXERR  defines  the  maximum 
deviation  in  the  turbulence  factors  between  iterations.  It  is  assumed  initially,  in  the 
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first  step  of  the  iterative  scheme,  that  the  flow  is  in  the  laminar  regime,  that  is,  all 
turbulmce  factors  =  1 . 

6.  Input/Output:  Parameters,  Conventions, 
Procedures 

It  follows  for  the  definition  of  specific  weight  that 

7  =  Pff  (46) 

where  g,  the  gravitational  constant,  is  32  ft/sec*.  Using  eq  (46)  the  computer  program 
computes  the  density  which  is  needed  for  calculation  of  the  Reynolds  numbers  eq  (8). 

Figure  5  defines  the  convention  for  distinguishing  the  entry  (F)  and  exit  (T)  nodes 
of  a  particular  branch.  The  arrow  defines  the  direction  of  positive  flow. 


F  T 

*  > 

Figure  5:  Conv«ition  defining  nodes. 

The  required  inputs  to  the  FORTRAN  program  code  flow. for  provided  in  Ap¬ 
pendix  A  are  displayed  in  Table  I. 

The  calculations  are  done  in  formal  English  uxuts  of  (ft.  Lb,  sec).  However,  the 
inputs  and  outputs  are,  when  appropriate,  given  in  conventional  units.  Conversion 
factors  for  these  cases  are  shown  in  Table  2. 
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l^ble  1:  Input  Variables 


Symbol 

Units 

i2,o 

§3 

typical  3500 

MAXERR 

§5 

— 

typical  0.05 

HLCF 

§2 

— 

typical  60.0 

Tlx 

§5 

— 

depends  on  network 

b 

§5 

— 

depends  on  network 

§2 

10~®  Ib-sec/ft® 

typical  HjO,  3.75 

W 

§6 

Ib/ft® 

typical  H2O,  64.0 

node  specification  — F,  —T 

§6 

— 

depends  on  pipe 

L 

§2 

ft 

depends  on  pipe 

d 

§2 

inches 

depends  on  pipe 

e 

§2 

mil 

depends  on  pipe 

source  specification  Q, 

§5 

gal/min 

depends  on  pump(s) 

Ikble  2:  Conversion  Factors 


Description 

Symbol 

Formal 

Units 

Conventional 

Units 

Conversion 

Factor 

Flow  rate 

Q 

ft^/sec 

gal/min 

450 

Pressure 

P 

Ib/ft* 

psi 

1/144 

Roughness 

e 

ft 

mil 

12,000 

In  order  to  facilitate  the  process  of  batch  input,  the  code  has  been  designed  to 
read  the  input  file  flowinp.dat  rather  than  require  the  user  to  laboriously  submit 
the  data  in  the  interactive  mode.  In  this  way  a  maximum  amount  of  flexibility  can 
be  incorporated  into  the  program  without  creating  a  tedious  data  entry  procedure  for 
the  user  of  the  code.  Thus,  the  ‘‘user  friendly”  goal  has  been  kept  in  mind. 


17 


The  outputs  from  the  code  flow. for,  which  are  listed  by  branch  number  (BR), 
are  displayed  in  Table  3. 

Where  appropriate  both  text  symbol  and  fortran  symbol  have  been  specified  in 
Table  4.  For  example,  PF  and  PT  refer  to  node  pressures  at  entry  and  exit,  respec¬ 
tively.  In  some  practical  cases  it  may  be  desirable  to  employ  the  equivalent  of  an 
“ideal”  constant  flow  rate  pump,  that  is,  qi,  =  q,k  in  eq  (24).  This  is  readily  handled 
by  forcing  yt  — »  0  or  equivalently  rk  — »  oo.  Because  of  the  head  loss  coupling  factor 
rule-of-thumb  introduced  at  the  end  of  Section  3,  the  most  convenient  way  to  force 
the  pump  source  to  behave  ideally  is  to  let  dk—^0.  In  practice,  because  of  numerical 
instabilities,  if  yk  —  0,  it  is  recommended  that  the  user  set 

d*  =  10”^  X  min{d<}  »  =  1,2, ...,6  butt^Ar.  (47) 

A  hardcopy  of  the  input/output  for  the  piping  networks  tested  is  provided  in 
Appendix  B.  A  more  detailed  description  of  these  cases  is  covered  in  the  next  section. 


Table  3:  Output  Variable 


Symbol 

Units 

BR 

§5 

Q 

§5 

gal /min 

PF 

§6 

10~®  psi  (mpsi) 

PT 

§6 

10~®  psi  (mpsi) 

HEAD  =  H 

§5 

10"®  psi  (mpsi) 

ERR 

see  RHS  eq  (44) 

— 

Re 

§2 

— 

TBF  =  V* 

§2 

— 

Y 

§5 

(gal/min)/mpsi 
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7.  Examples 

General  Observations 

Several  general  observations  can  be  easily  confirmed  from  the  data  (see  Appendix  B 
for  all  the  case  examples  discussed  in  this  section).  Specifically, 

PF  —  PT  =  HEAD  {pp  —  pr  —  h*) 

and 

Q  =  Q,  +  Y  HEAD  (9*  =  q^k  +  Vkhk) 

In  all  the  examples  the  sink  node  of  the  pump  was  arbitrarily  taken  as  the  datum 
node  for  convenience  of  interpretation.  The  node  pressures  under  columns  PF  and 
PT  are  as  expected,  strictly  positive.  In  all  three  examples  some  or  all  of  the  pipeline 
segments  were  operating  in  the  nonlaminar  regime.  This  is  evident  by  checking  the 
TBF  columns  of  each  example  in  Appendix  B.  The  associated  Reynolds  number  is 
also  provided.  The  branch  errors  were  calculated  according  to  eq  (45).  It  is  easily 
checked  that  the  MAXERR  parameter,  line  2  of  the  input  set,  put  a  ceiling  on  these 
branch  error  values. 

CASE  A  —  A  four  branch  pipeline  network  with  an  ideal  source  (see  Fig  6a  and  fib). 
This  case  illustrates  the  generation  of  an  ideal  flow  source.  Note,  as  discussed,  this 
condition  can  be  obtained  by  making  the  pipe  diameter  of  the  source  branch  relatively 
small.  The  associated  Reynolds  number  and  turbulence  factors  will  be  artificially  high 
and  may  not,  as  represented  in  the  line  1  output  for  Case  A,  be  within  the  defined 
format. 
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I  2 


(b) 

Figure  6:  a)  4  branch  with  ideal  source;  b)  4  branch  flow  graph  for  a). 


CASE  B  — -  A  six  branch  pipeline  network  with  nonideal  source  (set  Fig  4a  and  4b). 
Note,  as  seen  here,  the  effect  of  the  nonzero.a^mittance  of  node  1  is  to  impede  the  flow 
of  the  source  into  the  rest  of  the  pipdine  circuit.  Also,  the  source  is  required  to  be 
negative  because  the  pump  is  physically  driving  the  flow  opposite  to  the  convention 
defined  by  the  arrow  of  branch  #1. 
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CASE  .C  —  A  12  branch  piping  network  for  shoebox  cooling  rack  (see  Fig  7).  Note 
that  several  of  the  branches  have  negative  Q  values.  The  physical  interpretation  is 
that  these  cases  have  a  fluid  flow  opposite  to  the  convention  defined  on  Figure  7. 


Figure  7:  12  branch  shoebox  cooling  rack. 


8.  Future  Enhancements  in  the  Modeling 

Although  the  work  described  herein  extends  the  Ellis  (1993)  thesis  by  including  the 
nonlaminar  flow  regime,  it  still  should  be  considered  as  preliminary.  For  example: 

•  The  network  formalism  discussed  herein  should  be  compared  with  the  bench¬ 
mark  Hardy  Cross  method  [Lindebivg  J1987)]  as  well  as  other  more  recently 
proposed  schemes  [Carnahan  and  Christensen  (1972),  Bending  and  Hutchinson 
(1973),  Gay  and  Preece  (1975,  1977)]  in  order  to  assess  the  relative  computa¬ 
tional  efficiencies 

•  The  next  major  modeling  development  in  the  algorithm  would  be  the  addition 
of  the  capability  of  evaluating  the  node  to  node  heat  transfers. 
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•  Using  the  theory  for  compressible  fluid  flow  [Vennard  (1961)]  required,  for  ex- 
sunple,  when  the  cooling  fluid  is  a  gas,  the  developed  scheme  could  be  extended 
to  allow  for  compressible  fluids.  Due  to  the  sensitivity  of  gas  properties  to  tem¬ 
perature  this  development  would  have  to  include  effects  due  to  heat  transfer. 

•  An  altonate  but  similar  mathematical  formalism  could  be  derived  for  pumps 
best  characterized  as  constant  pressure  sources.  It  is  certainly  possible  to 
“Theveninize”  the  source  in  the  context  of  the  present  mathematical  formal¬ 
ism.  However,  it  should  be  noted  that  because  of  the  nonhnearity  inherent  in 
the  model,  this  Theveninization  would  need  to  be  repeated  iteratively,  decreas¬ 
ing  the  efficiency  of  the  process. 

•  As  previously  noted  in  Section  3,  a  significant  improvement  in  the  rate  of  con¬ 
vergence  of  the  method  was  obtained  by  reducing  the  discontinuity  in  the  fric¬ 
tion  factor  derivatives  at  the  transition  edges.  It  should  be  possible  to  apply 
mathematical  adjustments  on  eq  (19),  for  example,  using  a  technique  known  as 
Hermiteor  "osculating”  interpolation  [Hamming  (1962)],  in  order  to  completely 
eliminate  the  slope  discontinuity.  The  potential  gain  in  the  rate  of  convergence 
could  be  significant. 

•  Various  mundane,  but  no  doubt  practical  embellishments  could  also  be  consid¬ 
ered.  For  example,  an  option  of  noetric  system  units  could  be  added.  Also,  an 
option  to  specify  pipe  bends,  vertical  elevation,  and  cross-section  types  could  be 
added.  If  commercial  viability  is  of  concern,  the  program  could  be  dovetailed 
with  CAD  software. 

•  The  required  backgroimd  investigation  revealed  that  despite  a  long  history  of 
developments  on  modeling  fluid  flow  in  pipes,  only  recently  have  mathematical 
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methods  involving  the  theory  of  chaos  begun  to  unravel  and  predict  the  prop¬ 
erties  of  nonlaminar  fluid  flow.  A  small-scale  physical  model  which  predicts  all 
the  measured  properties  on  the  Stanton  diagram  is  at  this  date  unavailable. 

9.  Conclusion 

A  mathematical  model  for  predicting  laminar  and  nonlaminar  flows  in  pipeline  net¬ 
works  has  been  proposed  and  successfully  tested.  This  scheme,  upon  iterative  appli¬ 
cation,  has  been  found  to  converge  to  an  approximate  solution.  The  computational 
efficiency  of  this  convergence  process  was  found  to  be  highly  dependent  on  the  fric¬ 
tion  factor  modeling  in  the  transition  region  between  laminar  and  turbulent  behavior. 
In  particular,  a  significant  improvement  in  the  convergence  rate  was  obtained  by  re¬ 
placing  the  previously  proposed  linear  model  with  a  nonlinear  model  having  a  less 
dramatic  first  derivative  discontinuity.  A  linear  gauge  for  the  deviation  from  laminar 
behavior,  referred  to  as  the  turbulence  factor,  was  introduced  in  this  report  in  order 
to  facilitate  the  computational  task. 
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10:  Nomenclature 


Roman  Letter  Symbols 

A  pipe  cross-sectional  area 

d  pipe  diameter 

e  pipe  roughness 

/  friction  factor 
g  gravitational  acceleration 
h  head  loss  due  to  friction 
P  pressure 
Q  flow  rate 

R  pipe  resistance 

Re  Reynolds  number 
V  velocity 

Wt  time-rate-change  of  work  output  (pump) 
y  pipe  admittance 
z  elevation 

Greek  Letter  Symbols 

7  specific  weight 

e  pipe  roughness  ratio 

fi  viscosity 
p  density 


Units 

ft^ 

ft 

ft 

ft /sec* 

Ib/ft* 

Ib/ft* 

ft®/sec 

Lb-sec/ft® 

ft /sec 
Ib-ft/sec 
ft*/ (lb-sec) 
ft 


Ib/ft* 

lb-sec/ft* 

slug/ft* 
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Appendix  A:  Cholesky  Reduction 


This  appendix  provides  a  more  complete  description  on  the  Cholesky  reduction  algo¬ 
rithm  [Kraus  (1992)].  Positive  definite,  symmetric  matrices  can  be  factored  according 
to  the  rule 

M  =  LL^  (Al) 

where  L  is  an  upper  triangular  matrix.  Because  of  the  applicability  of  this  theorem 
[Ellis  (1993),  Kraus  (1987)]  to  the  network  matrix  Y„  eq  (41b)  it  follows  that 

Y„  =  LL^  (A2) 


and  for  eq  (40) 


LL^P  =  Q 


(A3) 


The  advantage  of  the  factorization  eq  (Al)  is  now  apparent  since  eq  (A3)  can  be 
broken  into  two  parts:  First, 

LJ  =  Q  (A4a) 


and  second, 

L^P  =  J  (A4b) 

which  are  easily  solved  in  succession  via  Gaussian  elimination  to  determine  first  the 
mliimn  vector  J  as  an  intermediate  step  then  the  derived  node  pressures  P. 
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Appendix  B;  Input/Output  Hardcopy 


****************  IIPUT  FORMAT  ♦*******♦••♦*•••*•♦•••*•• 

REYIOLOS’  lUNBER  FOR  RIGHT  EDGE  OF  TRIMSITIOH  REGIOH 

ERROR  TOLERAMCE  II  CALCULATIOH 

THE  HEAOLOSS  COUPLIIG  FACTOR  (TYPICAL  lUNBER  60  ) 

THE  HUMBER  OF  lODES 
THE  lUMBER  OF  BRAHCHES 
THE  VISCOSITY  Z  10  -6  LB-SEC/FT*3 
THE  DEISITY  LB/FT  “3 

STARTIHG  lODE.  EIDIHG  lOOE.  LEIGTH (FT) .DIAMETER (IH) ,  ROUGHHESS(MILS) 

lOTE  FOR  IDEAL  CURREHT  SOURCE  SET  DIAM-.OOIX  SMALLEST 
THE  lUMBER  OF  SOURCES 
SOURCE  STREMGTH  GAL/MIH 
SOURCE  BRAICH 


IIPUT  (CASE  A  4  BRAICHES  IDEAL  SOURCE) 
3S00.0 
0.005 
60.0 

4 

4 

3.76 

64. 

4.1.  1.0. .001.1.0 

1.2.  1.0.1. 0.2.0 

2.3.  1.0. 1.0.0. 5 

3.4.  1.0. 1.0.0. 6 
1 

6.0 

1 


OUTPUT  CASE  A 


BR 

Q 

PF 

PT 

BEAD 

ERR  R* 

TBF 

Y 

1 

6.00 

.0 

282.1 

-282.1 

»**«*4i 

.0000 

2 

6.00 

282.1 

183.8 

98.2 

.0032  10806.2 

5.6 

.0609 

3 

6.00 

183.8 

92.1 

91.7 

.0032  10806.2 

5.2 

.0652 

4 

6.00 

92.1 

.0 

92.1 

.0032  10806.2 

5.3 

.0649 

IIPUT  (CASE  B  6  BRAICH  PIPELIIE  lETVORK) 
3500.0 
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0.01 

60.0 

5 

6 

3.75 

64. 

1.5.  1.0. 1.0. 1.0 

1.2.  1.0.1. 0.2.0 

2.3.  1.0. 1.0.0. 5 

2.4.  1.0. 1.0, 0.6 

3.4.  1.0. 1.0. 0.6 

4.5.  1.0. 1.0. 5.0 
1 

-6.0 

1 

OUTPUT  CASE  B 


BR 

Q 

PF 

PT 

HEAD 

ERR 

R« 

TBF 

Y 

1 

-1.75 

25.1 

.0 

25.1 

.0073 

3149.7 

2.0 

.1679 

2 

1.75 

25.1 

14.6 

10.6 

.0075 

3149.7 

2.1 

.1645 

3 

.58 

14.6 

12.9 

1.7 

.0000 

1049.9 

1.0 

.3409 

4 

1.17 

14.6 

11.2 

3.4 

.0000 

2099.8 

1.0 

.3409 

5 

.58 

12.9 

11.2 

1.7 

.0000 

1049.9 

1.0 

.3409 

6 

1.75 

11.2 

.0 

11.2 

.0078 

3149.7 

2.2 

.1554 

*««««*«***************«**4ii»i*«**«**4i4i*****«>******4>*4i************4i4>***4i 

CASE  C  (  12  BRANCH  SHOEBOX  COOLING  RACK) 

3500.0 

0.01 

60.0 

8 

12 

3.75 

64. 

8.2.  1.0. 1.0. 1.0 

2.3.  1.0.1. 0.2.0 

3.4.  1.0. 1.0.0. 5 

4.8.  1.0. 1.0.0. 6 

5.8.  1.0. 1.0.0. 6 

2.6.  1.0.1. 0.5.0 

7.3.  1.0. 1.0. 1.0 

1.4.  1.0. 1.0.2. 0 

6.5.  1.0. 1.0.0. 5 

6.7.  1.0. 1.0.0. 6 
7.1.  1.0. 1.0.0. 6 
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5.1,  1.0, 1.0. 5.0 
1 

6.0 

1 

OUTPUT  CASE  C 


BR 

q 

PF 

PT 

HEAD 

ERR 

Re 

TBF 

Y 

1 

3.48 

.0 

24.6 

-24.6 

.0341 

6273.0 

3.4 

.0990 

2 

1.76 

24.6 

14.5 

10.1 

.0393 

3168.5 

2.0 

.1681 

3 

1.33 

14.5 

9.6 

4.9 

.0084 

2394.5 

1.3 

.2685 

4 

1.74 

9.6 

.0 

9.6 

.0371 

3138.2 

1.9 

.1750 

5 

1.74 

9.6 

.0 

9.6 

.0366 

3134.8 

1.9 

.1752 

6 

1.72 

24.6 

14.4 

10.2 

.0386 

3104.5 

2.1 

.1635 

7 

-.43 

13.2 

14.5 

-1.3 

.0000 

773.9 

1.0 

.3409 

8 

.41 

10.8 

9.6 

1.2 

.0000 

743.6 

1.0 

.3409 

9 

1.32 

14.4 

9.6 

4.8 

.0072 

2381.3 

1.3 

.2710 

10 

.40 

14.4 

13.2 

1.2 

.0000 

723.2 

1.0 

.3409 

11 

.83 

13.2 

10.8 

2.4 

.0000 

1497.1 

1.0 

.3409 

12 

-.42 

9.6 

10.8 

-1.2 

.0000 

753.5 

1.0 

.3409 

n  o 


Appendix  C:  Computer  Code 


Idabug 

ILIST 

c  prograa  for  CALCUUTIIG  PRESSURES  AID  FLOW  RATES  IM  PIPILIIE  lETWORKS 

C  a****************************************************************** 

C  PROGRAM  flov.FOR  WRITTEI  BY  PROF.  RON  J  PIEPER 


C  IPS  408  6562101 

C  AUG.  30.  1994 

C  SEE  FILE  flovIIP.DAT  FOR  IIPUT 

C  SEE  FILE  flovOUT.DAT  FOR  OUTPUT 

C  SEE  FILE  flovDIA.DAT  FOR  DIAGIOSTICS  01  IIPUT 

C  DESIGIED  TO  HAIDLE  UP  TO  40  BRAICHES 

C  CAI  BE  ADJUSTED  BY  IICREASIIG  DINEISIOI 

C  OF  THE  ARRAYS 

C  a******************************************************************* 

c  mu  fluid  viscosity  (lb-ssc/ft"2) 

c  rho  fluid  dansity  (lb/ft"3) 

c  g  accalaration  dua  to  gravity  (ft/sac‘2) 

c  nu  spacific  gravity  rhozg 

c  pi  pi 

c  rani  raynold’s  i  (rho*val*d/mu) 

c  f  friction  factor 

c  val  valocity  of  fluid 

c  pras  prassura  of  tha  fluid 

c  ****a*a*«***aa***pip«  variablas  ************************************** 
c  all  pipa  langth(s) 

c  d  pip*  diamatar(s) 

c  a  roughnass 

c  rat  roughnass  ratio 

C  HCLF  HEADLOSS  COUPLIIG  FACTOR 

c  **********  F  FACTOR  a******************************* 

C  TURB  TURBULEICE  FACTOR  -  F*2100/64  (-1.0  LAMIIAR  REG) 

C  NAXERR  NAXERR  H  TURBULEICE  FACTOR 

c  FOR  REFEREICE  ALSO  SEE  ITERATIVE  COOLBROOKE-HHITE 

C  AID  ALSO  SEE  DIRECT  CHURCHILL  aQ  (  II  PLOT  PROGRAM) 

c  spplias  for  rant  >  REYEDGE 

hedge  the  left  edge  of  TRAISmOl  REGIOI 

R^mrOLDSi  WHERE  COOLBROOE  BE6IIS 

c  ************************************************************** 
c  aaaaaa  natvork  constants  a*a*aaa****a*aa**********aa***a*****a 
c  I  tha  numvar  of  junctions  in  tha  systam 
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th«  •  junction  ~  datum  noda  (n-1) 
th«  i  branchas  in  tha  syatam 
tha  •  of  tourcaa 


c  II 
c  B 
C  S 


c 

c 

C 

C 

C 

c 

c 

C 


ASSUMPTIOIS 

1.  COOLBROOKE-WHITE  FORMULAE  FOR  REYi>  radga 

2.  HORK  II  EIGLISB  UIITS 

3.  FLUID  IS  nCOMPRESSIBLE 

4.  ALL  PIPES  ARE  CIRCULAR  II  SHAPE 

5.  THE  CORRECT  SOLUTIOI  IS  OBTAIMED  BY  ITERATIOH 


intagar  l.ll.B.QB.S.OI 
intagar  IT(40),  IF(40),X,  itast 

raal  mu.  rho,  g,  pi,  C(40.40).  d(40),  all(40),  r(40) ,QS1,QS(40,1) 
REAL  RLAM(40),TURB(40).Ct(40.40).Cy(40.40).YLAM(40),  Y(40.40) 
REAL  IS(40.1),YH(40.40),P(40.1),  YH(40,1) .V(40.1) ,  REY(40,1) 

REAL  AREA(40) .ISI(40.i) .YII(40.40) .Q(40.1) .ERR(40) .ap(40} .rat(40) 
REAL  TURBF.MAX£RR,qG(40).FCROS(40) 

REAL  REDGE.  EPS.H(40,i) .HCLF 
CHARACTER*20  FORM 
CHARACTER«3  BR 
CHARACTER*?  lAME(ll) 

BR-'BR  * 

HAME(7)«»  TBF  * 

IANE(4)«’  BEAD* 

IAI!E(2)-*  PF  * 

IAME(3)-*  PT  * 

IAME(6)-'  Ra  ’ 
lANE(i)-*  q* 

HAME(8)«’  Y  ’ 
lANECS)-'  ERR’ 


4  FORMAT  (A12} 

OPEH  (UIIT-7 .FILE-’flovDIA.DAT’ .STA7US>*IEH> .ACCESS- ’ SEqUEITIAL * ) 
OPEI  (UlIT-B.FILE-’flovIlP.DAT*,  STATUS- ’OLD’ .ACCESS- ’SEqUEITIAL*) 
OPEI  (UIIT-9.nLE-’flo«0UT.DAT',  STATUS- *IEH* .ACCESS- ’SEqUEITIAL’) 


g-32.174 
pl-3. 1416926 
READCB.B) 

8  FORMAT 

C  «rita(*,13) 

writa(7.13) 

13  format  (’IIPUT  THE  REYIOLDS  •  FOR  RIGHT  EDGE  TRAISITIOH  ’) 


30 


r«ad(8,30)  radge 
URITE(7.30)  HEDGE 
C  vrite(*,23) 

write (7, 23) 

23  fozvat  ('  enter  the  aaziauB  error  II  TURBULEICE  FACTOR  desired'  ) 

read(8.30)  NAIERR 
iniITE(7.30}  NAXERR 
C  WIITE(*,21) 

21  FORMAT  CEITER  THE  AVERAGE  HEADLOSS  COUPLIIG  FACTOR.TYPICAL  #60 ’} 

REA0(8.30)  HLCF 
write(7,21) 
write(7.30)  HLCF 
C  write (*.10) 

write(7,10) 

10  format  ('  enter  the  number  of  nodes  in  the  rack  system  l>100  ’  ) 
read(8.15)  I 

Il-M-1 

IS  format  (IS) 

write(7.1S)  I 
C  writa(*.lS)  I 

C  write (*,20) 

write (7, 20) 

20  format  ('  enter  the  number  of  branches  in  the  rack  system  B>H  *  ) 
read(8.1S)  B 
C  write(*,15)  B 

write(7,lS)  B 


C  write(*,2S) 

write(7.2S) 

2S  format  (’  input  the  wiscosity  _  X  10“-B  lb-SEC/ft“2  '  ) 

read(8,30)  mu 
VRITE(7.30)  mu 
C  tfRITE(*,30)  mu 

mn"mn*. 00001 
30  format  (flO.4) 

C  write(*,3S) 

write(7.3S) 

3S  format (  *  input  the  specific  weight  lb/ft*3*  ) 
read(8,30)  sw 
rhoasw/g 

C  write(*,30)  sw 

write(7,30)  sw 
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lode-branclx  aatriz 


c 
c 

c  utrlz  initialization 
c 

do  40  i«l.ll 
do  40  j-l.B 
40  C(i,j)-0 

do  42  i-l.B 
IF(i)-0 
42  ■T(i)-0 

C 
C 

C  racaiYa  dat  froB  tha  kaybord  to  davalop  a.d  and  all  Batricas  D 

c  for  circtilar  pasagas 

c  start  loop  on  branchas  to  input  noda  to  node  infonnation 
do  50  i-l.B 
*rita(7,55)  i 

55  forBat  (’  at  branch  nuBbar  *.  2z,  i5) 

C  Brita(*,60) 

HRLTE(7.60) 

60  foraatC* ENTER  nf(i).nt(i),langth(ft).diaBatar(in),ROUGHNESS-'Bil’) 

raadCS,*)  IF(i) ,IT(i) ,all(i) ,d(i) ,ap(i) 
ap(i)-ap(i)a.001 
rat(i)-ap(i)/d(i) 
c  convart  froB  Bils  to  inches 


WRITE(7,*)  IF(i)  .IT(i)  ,all(i)  .d(i)  ,^(i)  ,rat(i) 

C  WIITE(*,*)  IF(i),IT(i),all(i),d(i),ap(i),rat(i) 

C  70  fornatC  i5.i5,f8.3.f8.3,f9.5) 

c  *************  calculate  tha  coolbrook  crossing  value  at  R-  3500 
c  colabrook-vhita  formula  page  198  of  fox,  let  f  goto  4f 

sqf- (64. 0/2100.0} a*. 5 
3100  sava-sqf 


sqf-1.0/(-2.0*logl0(2.51/(REDGE*sava)  ♦  rat(i)/3.7)) 
IF  (SAVE  .EQ.  SQF)  GOTO  4000 
arror-abs (sava-sqf ) /sqf 
if  (error  .gt.  .001)  goto  3100 
4000  fcros(I)-sqf**2 

C  VRITE(*.«)  ’I.FCR0S(I).RAT(Z)'.  I.  FCROS(I).  RAT(I) 

50  continue 

C  SECTION  ON  SOURCE  INPUT 
c  INPUT  SOURCES.  S  OF  THEN 

C  qS  MATRIX  OF  FLOH  SOURCES 

C  QB  SOURCE  BRUCH 
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C  QSl  SOURCE  STRENGTH 

C  ilRITE(*,90) 

WRITE(7,90) 

90  FORMAT (  *  INPUT  THE  NUMBER  OF  SOURCES  *) 

REA0(8.92)  S 
writ* (7. 92)  S 
C  writ«(*,92}  S 

92  FORMATCIS) 

DO  98  I-1,S 
C  «IITE(*.95)  I 
WIITE(7,9S)  I 

95  FORMATC  'INPUT  THE  SOURCE  STRENGTH  FOR  THE’,  IS.'TH  SOURCE') 
READ(8,30)  QSl 
writ* (7. 30)  QSl 
C  vrit«(*.30)  QSl 

QS1-QS1/4S0 

C  CONVERT  FROM  GAL/MIN  TO  FT^3/SEC  CONVERSION  FACTOR  1/450 
C  X  GAL/MIR>  (450)--!  FT‘3/SEC 
.  WIITE(7,93) 

C  iniITE(*,93) 

93  FORMATC  INPUT  THE  BRANCH  NUMBER  OF  THE  SOURCE  ') 

REA0(8,15)  QB 
writ* (7. 15)  QB 
C  ¥rit«(«,15)  QB 

98  CONTINUE 

DO  100  K-1,B 
IF  (K  .EQ.  QB)  THEN 
QS(K.l)-  QSl 
ELSE 

QS(K,l)-0 
ENDIF 
100  CONTINUE 

c  Start  satting  up  ’  C  '  aatriz 
do  120  i*l,B 

if  (NF(i)  .LT.  N)  than 
C(NF(i),i)-l 

andif 

ifC  NT(i)  .LT.  N)  than 
C(IT(i),i)— 1 

andif 

c  conrart  diaaatar  in  inchas  to  faat 
d(i)«d(l)/12.0 
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o  rj 


«r«a(i)>pi*(d(i)**2)/4 

c  60  diaaatars  accounts  for  in/out  hsad  loss 
addl-HLCF*d(i) 

oll(i)*langth4  addl 

C  CALCULATE  THE  LAMIIAR  RESISTANCES  FOR  EACH  BRANCH 
RLAM(i)-  32.0*  ni«  *ll(i)/(AREA(I)  *  (d(i))**2) 
YLAM(I)«1/RLAM(I) 

120  CONTINUE 

C  HRITE  "C"  MATRIX  TO  SCREEN  AND  RECORD  FOR  CHECK 
imiTE(7,123) 

C  HRITE(*.123) 

123  FORMATC  *  THE  C  MATRIX  BELOV.  TRANSPOSED  '  ) 

DO  126  I-1,B 

C  VRITE(*,12S)  (C(J,I).J-1,N1) 

VRITE(7.125)  (C(J,I),J«1,H1) 

125  F0RMAT(8(1X.F8.4)) 

126  CONTINUE 

C  VRITE  THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS 

HRITE(7.*)  ’THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS’ 

C  HRITEC*,*)  ’THE  LAMINAR  RESISTANCE  AND  ADMITTANCE  VECTORS’ 

DO  130  I«1.B 

VRITE(7,*)  RLAM(I).YLAM(I) 

C  VRITEC*,*)  RLAM(I) ,YLAM(I) 

TURB(I)-1.0 
130  CONTINUE 

MATRIX  MULTIPLICATION  OF  MATRICES 
C 

DO  140  J-1,B 
DO  140  K-l.B 
Y(J,K)-C 

IF  (J  .EQ.  K)  THEN 
Y(J.J)  •  YLAM(J) 

EHDIF 

140  CONTINUE 

C  MATRIX  MANIPULATION  SECTION 
c  N1  rows  of  C 
c  b  coltams  of  C 
c  CT  transposo  of  C 

CALL  TRANSPCC.CT.Nl.B) 

C  writsC*.*)  ’  writs  transposs  of  watrix  C  * 


34 


writ* (7,*)  ’  writ*  transposs  of  matrix  C  ' 
call  WarrayCCT.B.Ml) 


C  mrita(*,*)  *  vrita  matrix  Y’ 
writa(7,*)  *  writa  matrix  Y* 
call  Harray(Y,B,B} 


C  PROPOSED  EKTRY  POIRT  FOR  ITERATIVE  SCHEME 
2S0  call  matmul(Cy,C.Y.Il,B,B) 

C  vritaC*.*}  *  vrita  PRODUCT  C*Y  IIXB  ' 
vrita(7,*)  *  vrita  PRODUCT  C*Y  IIXB  ’ 
call  Harray(Cy.Il.B) 

CALL  NATMUL(YN.Cy.Ct.ll.B,ll) 
c  noda  flov  sourca  vactor  CQs 

call  matauK  la.C.QS.Ml.B.l) 
c  commantad  out  artifact  of  old  achama 
DO  200  I-l.Rl 
ISd.D— IS(I,1) 

ISI(I,1)-IS(I,1) 

200  COITIMUE 

DO  29S  I>1.H1 
DO  295  J-1,11 

YII(I.J)»YI(I,J) 

295  COITIIUE 

C  SECTIOI  ALSO  WRITES  THE  ARRAY  PRIOR  TO  BEIHG  cHOLESKY  PROCESSED 
HRITE(7,*)  'YMI  MATRIX  PRIOR  TO  CHOLESKY  ' 

C  WRITEC*.*)  *YII  MATRIX  PRIOR  TO  CHOLESKY  * 
call  WarrayCYlI.ll.ll) 

WRITE (7, ♦)  »ISI  MATRIX  PRIOR  TO  CHOLESKY  » 

C  WRITEC*,*)  ’ISI  MATRIX  PRIOR  TO  CHOLESKY  ' 
call  WarrayCiSI.P* ,1) 
call  CHOLESKY(YII.ISI.Il) 

WRITE(7.*)  'YII  MATRIX  FOLLOWIIG  CHOLESKY  ' 

C  WRITEC*,*)  *TII  MATRIX  FOLLOWIIG  CHOLESKY  * 
call  WarrayCTII,ll,ll) 

C  P  lODE  MATRIX  P-YI“-1*IS 
C  H  BRAICB  PRESSURE  CT*P 
C  Q  BRAICH  FLOW  RATE 
C  TH  MATRIX  PRODUCT  T*R 
DO  296  1>1,11 
PCI.D-ISICI.I) 

296  COITIIUE 
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CALL  MATNULCH.CT.P.B.Ml.l) 

CAU  NATNUL(YH.Y.H.B.B.l) 

DO  300  I-l.B 
Q(I,l)-YH(I.l)+QS(I.l) 

C  DIVIDE  PRESSURES  IN  LB/FT‘2  BY  144  TO  CONVERT  TO  PSI 
300  CONTINUE 

c  REYNOLDS  NUMBER  REGIME  DETERMINATION 
DO  5000  I-l.B 

V(I,1)-ABS(  Q(I.1)/AREA(I)) 

REY(I . l)-rho*V(I . l)*d(I)/BU 
REN-REY(I.l) 

C  HRITEC*.*)  *  REYNOLDS  REY(I.l) 

WRITE(7,*)  *  REYNOLDS  t*.  REY(I,1) 
if  (REY(I.l)  .LE.  2100.0)  THEN 
YCI.D-  YLAM(I) 

TURB(I)-1.0 
TURBAV-l.O 
ERR(I)-0 
F-54. 0/2100.0 
•Isa 

if  (REYd.l)  .IE.  REDGE  }  THEN 
Xl-REDGE-2100.0 
Yl-Fcros(i)  -  64,0/2100.0 

F-  64.0/2100.0  ♦  Y1*SIH(PI*  (REN-2100) /(2*I1)) 

•ls« 

c  col«brook-vhit«  formula 

sqf- (64. 0/2100.0) •♦.6 
1310  savo-sqf 

sqf-1.0/(-2.0*logl0(2.Sl/(R«n*saT«)  *  rat(i)/3.7)) 

C  IF  (SAVE  .EQ.  SQF)  GOTO  1400 

•rror-abs (sar«-sqf ) /sqf 
if  (•rror  .gt.  .001)  goto  1310 
1400  f-sqf**2 

•ndif 

TURBF-F*R«n/64.0 

TURBAV- (TURB (I) ♦TURBF) /2 . 0 

Y(I,I)-1/(RLAM(I)*TURBAV) 

ERR(I)-ABS(TORBAV-TURB(I))/turb(i) 

TURB(I)-TORBAV 

ENDIF 

C  iniITE(*.*)  I. ’FRICTION  FACTOR’.  F.’THE  TURB  FACTOR’ .TURBAV 
imiTE(7.*)  I. ’FRICTION  FACTOR’,  F.’THE  TURB  FACTOR’ .TURBAV 
5000  CONTINUE 
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DO  5100  >1,B 

C  HRITEC*.*)  *  THE  ERROR  FOR  BRANCH '.J.*  IS*.ERR(J) 

HRITECT.*)  *  THE  ERROR  FOR  BRANCH*. J.’  IS*.ERR(J) 

IF  (ERR(J)  .GT.  EPS  )  EPS>ERR(J) 

5100  CONTINUE 

IF  (EPS  .GT.  HAXERR)  GOTO  250 

C  C  CCCCCCC  NOTE  ERR  DEFINED  IN  TERNS  OF  TURBF  IS  ALREADY  NORMALIZED 
VRITE(9,5900)  BR.NAME(l) .NAME(2) .NAME(3} .RANE(4) . 

*  NANE(5).NAME(6).RAME(7).NAME(8) 

5900  F0RMAT(A3.8(1X,A7)) 

C 

C  READY  TO  BEGIN  OUTPUT  OF  SOLUTION 

C  VRITE(*.*)  *BRANCH.FLOU(Q).PF(MPSI).PT(MPSI).  head.  ERR. 

C  *  REYt.  TURB  FACTOR.  FLOH(Q-QS).  * 

HRITE(7.*)  'BRANCH. FLOV(Q).  presURE  F.  pras  T.  head.  ERR. 

4  REYf.  TURB  FACTOR.  FLOH(Q-Qa)> 

DO  7000  J-l.B 

c  Sourcaa  calculated  in  terms  of  effective  pressure  diff  (  thevenin) 
C  NOTE  THE  FLOW  VALUES  ARE  IN  GAL/KIN 
C  CONVERT  FLOH  VALUES  TO  GAL/MIN 

C  DIVIDE  PRESSURES  IN  LB/FT‘2  BY  144  TO  CONVERT  TO  PSI 
C  ADMITTANCE  Y  CONVERTED  TO  (GAL/MIN) /nPSI 
c  *  1000  to  convert  to  apsi 

PF-P(NF(J),1)/144.0*10**3 
PT-P(NT(J) ,1)/144.0*10**3 
H(j.l)»H(j.l)/144.0*10**3 
Y(J. J)-Y(J.J)*450. 0*144. 0/10e*3 

qs(j.i)-45o.o*qs(j.i) 
q(j.i)-q(j.i)*450 
c  qG(j)-q(j.i)-qa(J.i) 

C  WlITE(*.6000)J.q(J.l) .PF,PT.H(j .1) .ERR(j) . 

C  4  REY(J.1),T0RB(J).Y(J.J) 

HRITE(7,6000) J.q(J, 1) .PF,PT.H(j , 1) ,ERR(j) , 

4  REY(J.1),TURB(J).Y(J,J) 

WlITE(9,6000)J.q(J.l),PF,PT,H(j.l).ERR(j), 

4  REY(J,1).TDRB(J).Y(J,J) 

6000  F0RMAT(I3,1X,F7.2.3(1X.F7.1).1X.F7.4,2(1X.F7.1) ,1X,F7.4) 

7000  CONTINUE 

HRITE(*,*)  »  *♦*♦♦**•*•*♦♦♦♦**♦**♦♦***•♦*********•♦♦♦* 

«IITE(*,*)  ’  ' 

HRITE(*.*)  *  SEE  FILE  flovOUT.DAT  FOR  OUTPUT  > 

HRITE(*.*)  *  SEE  FILE  flovDIA.DAT  FOR  CHECKING  INPUTS  > 
HRITE(*.*)  ’ 


TOITEC*,*)  *  ♦♦♦•••*••♦•***•*••*•*•*•*••**♦**********• 

STOP 

EID 

aubroutins  cholaakyCa.b.n) 

C  SOLVES  THE  PROBLEM  AX«B 
C  X  UIKIOWI  I  TUPLE 
C  A  OOUM  RAHX  2  MATRIX  OF  ORDER  H 
C  B  OOHH  X  TUPLE  (  UPOM  EHTRY  ) 

C  B-X  (  UPOH  EXIT) 
raal  b(40,l) 
raal  a(40,40) 

vrita(7,*)  *  in  cbolaaky  routine,  natriz  ordar  n  follovs  * 

call  VarrayCa.n.n) 

call  dacon^(a,n) 

b(l,l)-b(l.l)/a(i,l) 

do  10  i«2,n 

d»b(i,l) 

il-i-1 

do  5  1-1. il 

S  d-d-a(l,i)*b(l,l) 

10  b(i,l)-d/a(i,i) 

b(n,l)-b(n,l)/a(n,n) 

nl-n-1 

do  30  1-1, nl 

k-n-1 

kl-k<fl 

do  20  J-kl.n 

20  b(k,l)-b(k.l)-a(k,j)*b(j,l) 

30  b<k,l)-b(k.l)/a(k,k) 

return 
end 

subroutine  BatnniKc.a.b.n.B.l) 

C  SOLVES  THE  PROBLEM  C-A*B 

c  C  UIKHOVH  MATRIX  Izn 

c  A  nOHI  MATRIX  NzH 

c  b  now  MATRIX  IXb 

real  a(40.40).  b(40.40}.  c(40.40} 

C  n  •  coluans  of  a  and  c 

c  ■  •  roes  of  a  and  coluans  of  b 

cl  •  rovs  of  b  and  c 

do  25  i-l,n 

do  25  J-1.1 
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25 


c(i,j)-0 
do  25  k>l,m 

c(i,j)»  c(i,j)+  a(i,k)*b(k,j) 

roturn 
•nd 

subroutin*  doconpCa.n) 
raal  a(40.40) 
i“l 

writa(7,*)  *  now  in  dacomp  routine  ' 
vrita(7,*)  *tha  order  of  aatriz  a  is*,  n 
do  410  ii-l.n 

write(7.35)  (a(ii. j j) ,jj-l,n) 

410  continue 

if  (  aCl.l)  )  1.  1.  3 

1  write(7.2) 

»rite(7,*)  i.aCi.i) 

2  fonatC  *  zero  or  negative  radicand  ’) 
goto  200 

3  a(l,l)>sqrt(a(l,i)) 
do  10  j>2,n 

10  a(l,j)»a(l,j)/a(l,l) 

do  40  i"2,n 
il«i-l 
d^adti) 
do  20  1-1. il 
dhold-d 

20  d-d-a(l,i)*a(l,i) 

if  (  d  .eq.  0}  then 

WRITE(7,e)  »  d  is  zero  ’ 
¥rite(7,*)  i.ad.i),  dhold 

endif 

srite(7,*)  ’  i,  aatrizCi.i)  follow  ’ 
wrlte(7,*)  i,  a(i,i) 
if  (aCi.i))  1,1,21 

21  a(i,i)-sqrt(d) 

if  (  i  .eq.  n)  goto  45 
i2-i+l 

do  40  J-i2,n 
d-a(i,J) 
do  30  1-1, il 
30  d-d-a(l,i)*a(l,j) 

40  a(i.J)-d/aCi.i) 

45  do  50  i-2,n 
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c  zero  out  lo««r  part  of  natriz 

do  50 

50  a(i.j)-0 

writ«(7,*)  *Bada  it  through  daconp  ,  matrix  follous  * 
do  400  iiBl.n 
vrita(7,35)  (a(ii,jj) 
format (8 (2z . f 8 . 4) ) 

400  continua 
200  ratum 
and 

SUBROtmiE  TRAISPU.AT.IR.fC) 

REAL  A(40.40).  AT(40.40} 

DO  20  >1.RC 
00  20  K-l.HR 
20  AT(J.K)-A(K.J) 

RETURH 

EHD 

subroutina  varrayCa.n.m) 

c  vritas  array  to  acraan  and  to  diagnostic  fila 
raal  a(40,40) 
c  n  rovs  of  a 

c  m  columns  of  a 

do  20  j«l,n 

writa(7.30)  (a(j ,i) ,i-l,m) 

C  arita(*,30)  (a(j ,i) ,i»l,m) 

20  continua 
30  format  (8(lz,f8.4)) 

ratum 
and 
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